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ABSTRACT 


A combined theoretical and empirical numerical model was 
developed which predicts the performance of continuous 
electrode coaxial magnetoplasmadynamic thrusters ав а 
function of thruster dimensions, mass flow rate, and input 
current. This model was used to predict the effects of 
scaling on these thrusters. 

The model predicts that for scaling factors down to one- 
half, relations can be found relating the performance of one 
thruster to another. The model was used to examine these 
relationships for four different thruster configurations 
over a broad range of operating currents. The thrusters 
examined consisted of two geometries and their half scale 
counterparts. A conclusion from the analysis is that 
scaling down the size of the thruster by 50% can reduce the 
total power input by 30% to 40% at comparable efficiencies. 
However, this is at the cost of increasing the specific 
impulse by a factor of two which may render the thruster 


inappropriate for the intended missions. 
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I. INTRODUCTION 


Equation 1 expresses the relationship between initial апа 
final space vehicle mass to the propellant exhaust velocity 
and the velocity increment required for a maneuver in space 
Ri... This equation shows that a rocket engine's exhaust 
velocity must be greater than the velocity increment of the 
mission of interest for a significant amount of payload to 
be delivered. If the exhaust velocity is much less than the 
velocity increment, the ratio of spacecraft masses becomes 
small implying a small final mass or payload. The velocity 
increment depends only on flight path and is independent of 
the vehicle. Therefore, the only way of improving the mass 
ratio is to increase the exhaust velocity of the rocket 


engine. 


Me Mi = exp(-Av/u,) (1) 


Chemical rocket engines have a maximum exhaust velocity 
determined by the amount of energy stored in the chemical 
bonds of the fuel, This limits exhaust velocities to less 
than 5000 meters per second. Orbit transfer missions such 
ПЕ оу Earth Orbit (LEO) to Geosynchronous Orbit (GEO) 
require velocity increments of approximately 4200 meters per 
second. (риса топ о: these values to Equation l gives a 
mass ratio of less than one-half. This implies that, at 


best, less than half of the initial weight of a spacecraft 


can be useful payload. For current space motors, the ratas 
of initial mass to final mass is about 4:1. An alterna е 
to chemical propulsion is electric propulsion. ЕЛес Ете 
propulsion avoids the thermodynamic limitations of chemical 
engines by employing an external power source to generate 
the forces required to accelerate a propellant. The amount 
of energy available to provide thrust is no Тоны 
determined by the nature of the propellants. This enables 
greater amounts of energy to be applied to the propellant 
than possible with chemicals. One such device is the 
magnetoplasmadynamic (MPD) thruster. This coaxial electric 
propulsion device uses the Lorentz jxB force to accelerate a 
plasma to velocities in excess of 25,000 meters per second. 
This increased exhaust velocity enables a 50% to 200% 
increase in the mass of the delivered payload. An electrical 
power supply must now be included in the spacecraft mass, 
but the reduction in propellant mass enabled by the MPD 
thruster's high exhaust velocities may be suffucient to 
allow for the mass of the power supply. Figure 1 shows an 


MPD thruster. 





Figure 1 Continuous Electrode MPD Thruster 


It can be seen that the thruster is structurally a very 
simple device. The MPD thruster has attracted considerable 
attention not only for its simplicity, but also because much 
greater thrust levels than any other electric propulsion 
device can be generated. This, combined with its ability to 
operate with virtually any propellant, make it a prime 
candidate for supporting space missions. 

The price to be paid for eliminating the dependence on 
ene chemical energy in the fuel and oxidizer is the 
requirement of an external power source. Power requirements 
for even modest sized MPD thrusters reach Lo the 


multimegawatt range. 


Space based power supplies of this size will not be 
available for many years. Two options are available to 
reduce the power requirements of the thruster. The first is 
to increase the efficiency of the thruster design. Since 
these thrusters are now operating in the 50% efficiency 
regime [2], this technique promises a reduction in power 
requirements of less than a factor of two. The second 
option is to develop smaller thrusters which require less 
power than their full size counterparts. MPD thrusters E 
virtue of their high thrust density, would still be able to 
produce useful levels of thrust after significant sca imen 
This technique promises much greater factors of power 
reduction than the first option because power requirements 
are related to size by a term of greater than linear scale. 

It is the goal of this thesis to develop a computer model 
of an MPD thruster which will predict specific impulse and 
efficiency as functions of current, thruster dimensions, and 


propellant mass flow rate. 


ШИШЕК ASNO TAE MPD THRUST MODEL 


Equations 2, 3, and 4 are the equations of conservation 


of mass, momentum, and energy respectively. 


9p/dt + У•(рч) = 0 (2) 


o[au/dt + (u-V)ü] = -Vp + (jxB) + fy (3) 


ШО ох т ое) (с т + 17/2) = dp/dt + jrE + Og + Oy - by (4) 


Since an MPD thruster is intended to be operated in a 
continuous mode, these equations can be simplified to the 
EN state forms of Equations 5, 6, and 7. This is 
justified because measurements of the start up and shut down 
transients in MPD thrusters are measured in seconds where as 


Киша! thruster operation is intended to last for up to 100 


days [3]. 


У•(рч) = 0 (055 


пе) = =ур = СјЈЕВУ t f. (6) 


o(uü*V)(cp5T + u*/2) = jeE + dg + фу – Qr vA 


For the девшуа то ог ап expression for the thrust ofan 
MPD thruster Equation 6, conservation of momentum, will be 
of paramount importance. The low density and viscosity of 
the gaseous argon, which will be the propellant used, within 
the arc chamber suggests that the viscous drag density can 
be assumed negligible. Calculations [4] show that Деши 
force, when integrated over a typical thruster, is only on 
the order of a tenth of a Newton where as the total thrust 
is typically three orders of magnitude greater than this. 


This simplification results in Equation 8. 


p(ü-V)ü = -Vp + ( jXB) (8) 


It is assumed the thruster in question is a small section 
of ап infinitely long coaxial device. This азосцши и 
greatly simplifies the mathematics by eliminating the 
requirement to deal with electric and magnetic field 
fringing and the exit plane of the thruster. In addition ШЕ 
will be assumed that the thruster consists of a series of 


coaxial streamtubes of cylindrical cross section and that 





steady state migration of propellant across the stream tube 
паге does not occur. This permits the application of 
one dimensional analysis to the problem of modeling the 
forces generated by the thruster. As will be seen later, 
this division of the thruster into coaxial streamtubes is 
also necessary for the application” of a numerical 
integration program, 
The above assumptions and simplifications permit Equation 

8 to be expanded using cylindrical coordinates to the form 
Equation 9. 

pl(ür+üg+u,)(rd /dr+69 /86+29 /dz)]-[ü,+ügtu,] = 

-[£8p/3r«683p/86«29p/8z]*[Cjr*3a* 32) (Bet Bo*Bz)] (9) 


From symmetry, the assumption that the only velocity 


EEUU Сап be maintained is in the axial or "z" direction, 


the definition of the vector dot product, Equation 9, 
eee pe Simplified to Equation 10. 


2002002792 - -fóp/Or - zop/Oóz + Zj,Bg – Рј„Вф (10) 


By equating vector components, Equation 10 can be split 


што Equations 11 and 12. 


pu ou dz - -ор/42 % 1186 (ЕЕ?) 


0 --Әр/дг - 1.8 (12) 


After rearranging, Equation 11 сап be put intor ni E K 


of Equation 13 and Equation 12 becomes Equation 14. 


f,n = u,au,/dz + dp/dz = jJ (Hop 


fp = -9p/or = j, Bg (14) 


It is convention [1], [5], І6|), [7], 18]; (9| (о ее 
to the f, component as the blowing force density and the I 
component as the pumping force density. The thrust produced 
by an MPD thruster is, in large part, the sum of these two 
components integrated over the volume of the thruster. The 
remaining thrust is the result of electrothermal heating moi 
the propellant as it flows through the arc chamber. 
Expressions for these components of thrust will be developed 
individually after which the three will be summed to give 
the total thrust. 

An mathematical model for thrust can only be accomplished 
if the current density within the arc chamber can be 
mathematically defined. ТЕ Һаѕ been experimentally 
determined [6] that over a wide range of thruster shapes and 
sizes, the current density follows an "S" shaped path as 


illustrated in Figure 2. 





ООО ТТГ 
CATHODE 


wile 2) Current Flow Pattern Inside an MPD Thruster 


For a continuous electrode thruster, the depth of this 
curve or amount of axial displacement is determined by the 
Hall parameter. For typical operating regimes, the Hall 
parameter varies from effectively zero to a value of 
approximately ten [5], [6] depending on the current for a 
Mesa thruster. By trial and error, Equations 15 to 18 
were developed and provide the current flow pattern 


required. 


for н п тег )/2 
Ss (1 37r]2- Jo 5225 10°(т-г)л°]*)° (15) 


E О па (16) 
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For (r 4ır,)/ ces 


jr (Or O СЕ (ІТ! 


A HO (18) 


With these relations to define the current within the 
thruster, expressions can be found for the blowing, pumping, 


and electrothermal components of thrust. 


BLOWING THRUST 


The blowing component of Br results from the vector 
product of radial current density with an azimuthalm 
oriented magnetic field as shown in Figure 1. Equations 15 
and 17 provide the radial current densities. The magnetic 
field is found by applying Ampere's law for straight 


conductors. 


В = 9В - $u$J/2mr Ampere's Law (19a) 


The magnetic field which results is given in equation (19b). 


10 





БК Ју спг)( 1->/2') (19b) 


As can be seen from an examination of this equation, the 
magnetic field is a function both of radial and axial 
position within the thruster. Vector multiplication of these 
two crossed vectors results in a vector directed in the 


positive "z" direction with a magnitude defined by Equations 


and 21. 


ИШ <г<(г+г„)/2 


ЖЕТ 12911,/2тс277(17/г17-(2.22х10%(г-г.)7212Ү7(1-2/2!)(20) 


Pee (rc+ra)/2<r<ra 


ШИЕ 5 -1597г,/2тг2!)(|7/г|7-12.22х10%(г,-г)27212)5(1-2/219(21) 


These last two equation- give the force density, in newtons 
Ee ubic meter, as a function of radial and axial position. 
To find the total blowing thrust component, it is necessary 
to integrate these force densities over the internal volume 


Sethe arc chamber. 


Eu ro«rC(rctra4)/2 


2 ‚(Fetra)/2 
Е,] = | \fzıdz2nrdr 220) 
0 


rc 
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For (тоатга) оа 


ка 
no E |f22d22"rar ( 233 
0 аа 


After integrating with respect to the axial and azimuti gm 
variables, Equations 24 and 25 result. 
а ЕЈ f 
Popa? - 95510" (elie - [2.2 jo та Ја н (24) 


Tc 


Га 


F2" |? 95x10* ([J/r]?-[2.22x109 (zg-1)J]* )5Jrdr (25 
Gen D 

The form of this equation is such as to resist integration 

analytically or through the use of generic forms given in 

integral tables. The total blowing force was obtained by 

applying numerical integration techniques to Equations 24 

and 25. This integration uses the program in Appendix B 


titled Theoretical Model Program. 


PUMPING COMPONENT 


The pumping force density given by Equation 14 is the 


12 


result of the vector product of the azimuthal magnetic field 
 Й ап axial component of the current in the negative "z" 
direction. This cross product gives rise to a force density 
component in the negative radial direction. This force 
density is balanced by a pressure gradient in the positive 
radial direction. Integration of the pressure gradient at 
ІНЕ 2 - 0 position yields an expression for pressure as a 
ШшшсЕетоп о: radial position. Multiplying pressure with a 
differential area and integrating provides the total pumping 
morce . 

Expanding Eequation 14 by substituting in the definitions 
of the axial current density and the azimuthal magnetic 
field at the z'- O0 location results in Equations 26 and 27. 
For ПА пе "ге )/2 


f£, --8p/8r = cJ*(r-r,)/z'r (26) 


For E rr.)/2&r&rg 


1E opor CUT —r)/2 E ( 295 


where c 2 7.08x10!" newtons/cm^amperes? 
separating variables and integrating gives Equations 28 and 
29 which collectively represent the equation for the radial 


pressure gradient as a function of radial position. 


ШОТ т <т<(г+г„)/2 
D po+c(r=r ¿1nr)J*/z' (0299 


13 


For (r tra) so na 


p = ПЕ ОЕ уя" (29) 


Integrating Equations 28 and 29, with respect to a 
differential ring area, and collecting terms yields equation 


30% 


Eoumpiag One 2 Sec een 
Ти (Бака) 2 ЈЕ точку оуу ыы = 
A A A O ns 
га ага ISO pr OI 


+re)/2)-ra((ratr.)/2)2/a-((r.+r.)/2)2/3)) 8 


while this expression is long, it is algebraically simple 
and is easily evaluated numerically. Like the blowing thrust 
term, Equation 30 has been included in the Appendix B 


progran. 
ELECTROTHERMA THRUST 


Due to ohmic heating the propellant in an MPD arc chamber 


increases in temperature from an initial 285 Ес 


14 


temperatures exceeding 20,000 K at some locations within the 
thruster. Cory [4] derived an expression EOE the 
electrothermal thrust based on the assumptions of negligible 
viscous forces, uniform pressure on the chamber walls, and 
no thruster pressure present on the outside surfaces of the 
arc chamber. His result for the pressure is presented here 


as Equation 31. 


р - 3.02х103,55%3 КӨЗЕ) 


Multiplying this expression by the exit plane area of the 
thruster provides an expression for the electrothermal 


component of thrust. This expression is Equation 32. 


ЕИО SM CPZ) 


ШООЛА THRUST 


ПСИ с ба Беа before, the total thrust is the sum of the 
hree thrust components. A Fortran program titled 
Theoretical Model Program given in Appendix B performs the 
Summation. The results of calculations using this model are 


sec ussed in section III. 


ks 


THRUSTER VOLTAGE 


An emperical expression for thruster voltage was 
assembled using all available information on coaxial MPD 
thrusters tested to date, [2], (3), [4], [5], [6]. [8] они 
and ll First, an equation for voltage as a function of 
thruster length was developed. Second, the influence of 
thruster diameter was included. Finally, the influence of 
mass flow rate was incorporated into the equation. This 
voltage model was tested using thruster geometries not used 
iN TES formulation. The model predicted interpolated 
voltages to within three percent and extrapolated voltages 
within seven percent. Both of these errors are within the 
plus-or-minus 10 to 30 percent error bars typical of 
laboratory testing associated with MPD thruster operation. 
Therefore, the voltage model was judged to be acceptable. 


This model is represented by Equation 33. 


V=((55.4-26.3]* (1-7, / 5a LA EZ IO) 
((39. 226558 *(1=r 75. DDr. Ele Ло E 
J+((35.7-54.5*(1-r3/5.1)+3.88*(1-2/21.6))*1E-8)*J2)* 


(178.0-47.3*1000.0*W+4.3*(1000.0*W)2)/49.0 (38) 


In this equation, z represents the thruster length and Wii 


propellant mass flow rate. 


SUMMARY 


From the input parameters of thruster length (z), anode 
radius un cathode radius «205 current (J), and mass flow 
rate (W), the model predicts thrust and voltage. These 
predictions are used in Equations 34 and 35 to calculate 


specific impulse (Isp) and efficiency (N) respectively. 


Шет - T/(W#Q B) (34) 


N = (0.5*W*(ISP*9,8)2/(V*J) (35) 


The model has been written as a Fortran program titled 


Theoretical Model Program. л серу 1S included, unir this 


thesis as Appendix B. 
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ІІІ. APPLICATION AND RESULTS OF THE THEGRET ECAR MOI 


The theoretical computer model was tested over a wide 
range of input currents of known thrusters. Graphs of 
efficiency versus specific impulse matched the experimental 
results [6] to within plus or minus ten percent. Similar 
agreement was obtained for graphs of efficiency versus input 
power. The plus or minus ten percent criteria was deemed 
acceptable because it placed the model's plots within the 
error bar of the experimental data. This test verified the 
model's accuracy. 

The model was used to predict the'specific impulse and 
efficiency of four different thrusters over a range of 
current levels. The thrusters were chosen with geometries 


which grouped them into two pairs of full and half sized 


devices, All thrusters were assumed to be using argon as 
propellant. A coded designator was created for each device 
to make identification easier. This designator is of the 


form A/B O/D: A is the thruster length. B is the anode 
radius. C is the cathode diameter. D is the mass flow rate 
of the argon propellant in kilograms per second. All 
dimensions are in centimeters. Thus a 20/5/17.006 thrustes 
is 20 centimeters long with a 5 centimeter inside anode 
radius and a 1 centimeter diameter cathode. The mass flow 


rate is 6 grams per second. Thruster sizes, aspect ratios, 
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mass flow rates, and applied currents are all consistent 
with those regimes of operation typical of MPD devices. 
Very small devices operating at power levels or mass flow 
rates where previous theoretical or experimental evidence 
indicated operation would be unstable were not examined. 

Ше Tour thrusters examined were 20/5/1/.006, 10/5/1/.006 
their half scale analogs 10/2.5/.5/.0015, 5/2.5/.5/.0015 
respectively. The 195/1/.006 configuration was chosen 
because it closely models the Princeton Benchmark thruster 
Ber which significant theoretical and experimental results 
exist. The 20/5/1/.006 thruster was selected because recent 
research indicates thrusters with extended electrodes are 
more efficient than shorter electrode thrusters like the 
Benchmark. The mass flow rate for the half scale thrusters 
was reduced by the square of |. factor to reflect the 
reduction in exhaust area. This maintains a constant atomic 
particle density. Failure to do so would place the scaled 
Mmimmusters in a completely different operating regime. 
Parameters for each of these devices became the input to the 
computer model with a range of currents from 1,000 to 25,000 
amperes. The model's outputs were plotted and are presented 
in Figures 3 through 7 inclusive. 

Figures 3 through 5 show the relationship between 


specific impulse and efficiency. A comparison between the 
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geraphs for 10/5/17.006and 2075/1792005, РОСО ССС usc 
that; as predicted, the longer configuration is more 
efficient at a given specific impulse. The same result is 
observed in Figure 4 for the half scale thrusters. This can 
be explained by examining the change in current density and 


current patterns as thruster lengths increase. For a given 


current, a longer device has a lower average current 
density. This results in lower collision frequencies and, 
therefore, smaller Hall currents. Hall currents, in MPD 


thrusters, travel in the axial direction and are the result 
of the combined effects of differing collision rates and 
mobilities between electrons and ions. It is the Hall 
current which creates the pumping component of thrust. 
Equation 18 is an emperically derived expression for the 
Hall current. Because this component involves propellant 
being accelerated through radia pressure gradients and a 90 
degree turn to exit the thruster, it is less efficient than 
the blowing component which is axially oriented throughout 
the acceleration process. This result could have been 
predicted by examining the relative magnitudes of the forces 
from equations О Шапа ро Qualitatively, both the 
blowing and pumping thrust components are functions of total 
current raised to the third power. However, the 
multiplicative constant for the blowing force appears to be 
larger than any value expected from the corresponding 
constant in the pumping equation. One would expect this 


phenomena to continue until viscous losses increase up to 
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large enough values to offset any gains realized by further 
increases in thruster length. These viscous loss effects 
were introduced to the model through the use of emperical 
data used in the formulation of the expression for thruster 
voltage. The dominance of the blowing component of thrust 
can also be observed by comparing the curves for thrusters 
10/5/1/.006 and ВО 25009570. 00 1588 Figure 5. These 
thrusters are the same length but the 10/2.5/.5/.0015 
Mister one half the radius of the 10/5/1/.006 thruster. 
This results, again, in an increased average current density 
and the corresponding decrease in efficiency for the smaller 
diameter thruster for constant specific impulse. 

Comparing the full size thrusters to their scaled 
counterparts discloses the scaled thrusters appear more 
efficient at a given specific impulse. This would seem to 
 шШгегайісес the explanation of current density losses in the 
previous paragraphs. However, examination of Figures 6 and 
7 show that the efficiencies of subscale thrusters. are 
MUSUDatic in nature as "opposed to more linear in the full 
EE UM hrusters when plotted against power. A continuation 
of the plots for the scaled thrusters to higher powers 
results in crossing of the full size and subscale thruster 
plots making subscale thrusters less efficient than the full 
sized thrusters. This behavior is indicative of the 
presence of a higher order loss mechanism beginning to 
dominate thruster operation at high specific power levels. 


MIS hisher Order term is most likely a reflection of frozen 


2: 


flow losses being introduced with the empirical elements 
used in the formulation of the computer model. Frozen flow 
losses are the result of energy stored in the excited states 
of the propellant as it is exhausted. As the propellant is 
heated and  ionized, energy is stored it the excited states. 
The propellant is exhausted while in an excited state 
resulting in an energy loss. The density of MPD thruster 
exhausts are too low to have collision rates high enough to 
de-energize these excited states and recoupe this energy in 
nozzles of any practical length. Frozen flow losses are 
related to the energy density of tthe propellant. Smaller 
thrusters operating at approximately ehe same power have 
greater amounts of energy per unit propellant mass and unit 
time because they have less mass іп which to deposit the 
energy. This effect is present throughout the range of 
operating conditions but is more pronounced at high power 
due to the increased occurrence of secondary ionizations. 
Figures 6 and 7 are plots of efficiency versus input 
power for the four thrusters examined. These graphs show 
that scaled thrusters can be expected to operate at lower 
power levels than their full sized counterparts along lines 
of constant efficiency. This phenomena, as the above 
explanation presents, is restricted to the lower power 
regimes of operation. The graphs show that overall power 
requirements can be reduced by using subscale thrusters. 
What is also shown is that for scaled thrusters operating 


along lines of constant efficiency, the smaller thrusters 
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have much higher specific impulses. This is not necessarily 
an additional benefit. Operationally, a propulsion system 
has a minimum total mass for a given total impulse 
requirement determined in part by the specific impulse [1]. 
Operating at specific impulses above the optimum point can 
result in an increase in propulsion system mass. This is 
the consequence of power supply mass exceeding the optimum 
value. The graphs show a one half reduction in linear scale 
results in a 30% to 40% reduction in power requirements for 
thruster operation in the 10% to 20% efficiency regime. 
This is accompanied by a 40% to 60% increase in specific 


impulse over the same range of efficiencies. 
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vee CONGCLUSTONS AND RECOMMENDATIONS FOR FUTURE RESEARCH 


The data generated by the computer model indicates scaled 
thrusters can operate at comparable efficiencies to full 
size thrusters but at reduced power levels. The reduction 
in power is on the order of 30% to 40% for linear 
geometrical scale factors of one half. While this reduces 
the power requirements for an operational system thereby 
making it more plausible, this reduction in power is 
accompanied by an even greater increase in specific impulse. 
The increase in specific impulse is on the order of 40% to 
60% putting the subscale MPD thruster into the 2500 sec 
regime. Studies [8] indicate orbit transfer missions are 
optimized with specific impulses in the 1600 second to 2000 
second regime. Scaling ос MPD Dans may result in 
operation at specific impulses higher than optimum, thereby 
reducing the amount of payload such a propulsion systen 
could deliver. This, coupled with the fact that this 
phenomena only occurs in the undesirably low efficiency 
range of 10% to 20%, indicates scaling thrusters to reduce 
power requirements may not be a good choice when developing 
an operational systen. 

Much work remains to be completed to verify this model. 
This will require extensive testing of a wide variety of 


different MPD thrusters and their scaled counterparts. An 
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attempt to do this is now underway in the Air Force 
Astronautics Laboratory's Electric Propulsion Laboratory. 
However, approximately two years of work remain before 
sufficient empirical data is generated to verify and further 
refine and expand the model. Once that is done, the model 
will be a valuable tool for the prediction of thruster 
performance. Appendix A is a data reduction program written 
in Fortan and intended for use in this Laboratory once the 
test apparatus is functional. An example of its output 


follows the program. 
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fat l1 “JN Ls a dd si 


DATA REDUCTION PROGRAM 


PROGRAM MPD 
CHARACTER HH*9,11*4 
БЕККИ (CAC) Оу /УА(20) ,МА(20),МА(20),51(407,5ТА 
1(20),Е(40),ҒА(20),Е(40),ЕА(20) 
WRITE(*,106) 
106 FORMAT(1x,'THIS PROGRAM ACCEPTS UP TO 33 DATA POINTS') 
984 CONTINUE 
WRITE(*,101) 
101 FORMAT(1X,'ENTER THE CATHODE LENGTH, ANODE LENGTH, ' 
E Tl CM) AND SCALE FACTOR IN F4.1') 
HESDC*,102» XX,YY,22 

102 FORMAT(3F4.1) 
WRITE(*,001) 

001 FORMAT(1X,'ENTER THE DATE AND TIME ACCORDING TO ' 

1,'AIR FORCE CONVENTION') 
ЕБЕШСЕ 002) НН,ТІ 
002 FORMAT(A9,A4) 
WRITE(*,981) 
981 FORMAT(IX,'ENTER THE MASS FLOW RATE IN KG/S IN F6.4') 
READ(*,086) DM - 
086 FORMAT(F6.4) 
DMM=DM*1000. 
WRITE(*,199) 

199 FORMAT(1X,'HIT CONTROL-P THEN RETURN TO PRINT INTRO.' 
1,'AFTER PRINTING, REPEAT THE COMMANDS AND REWIND THE' 
ШІ БІРЕК Т0 ТНЕ ТОР.") 

READ(*,198) КР 

198 FORMAT(A1) 

WRITE(*,200) 
200 FORMAT(9X,' b) 
BURITE(*, 201) 
201 FORMAT(9X,'| THIS DATA WAS OBTAINED FROM |") 
WRITE(*,202) 

202 FORMAT(9X,'| THE VARIABLE GEOMETRY MPD |') 
WRITE(*,203) 

203 FORMAT(9X,'| THRUSTER AT THE AIR FORCE |') 
WRITE(*,204) 

204 FORMAT(9X,'| ROCKET PROPULSION LABORATORY |') 
WRITE(*,003) HH,II 

003 EORMAT(9X,'!| ON ',A9,' AT ',A4,' HOURS. |') 
WRITE(*,205) 

205 FORMAT(9X,'| VOLTAGE (V) AND CURRENT (J) |') 
WRITE(*,206) 

206 FORMAT(9X,'| MEASUREMENTS WERE READ FROM AN |') 
WRITE(*,207) 

207 FORMAT(9X,'| OSCILLOSCOPE ATTACHED TO THE |") 
WRITE(*,208) 

208 FORMAT(9X,'| THRUSTER ELECTRODES. THE TOTAL |') 
WRITE(*,209) 
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209 
210 


211 


215 
214 
215 
216 
21 
218 
215 
2209 
221 
222 
223 
224 
оо 
226 
227 
228 
2 
230 
291 
400 


234 


FORMAT(9X, 
WRITE(*,21 
FORMAT(9X, 
WRITE(*,21 
FORMAT(9X, 
WRITE(*,21 


x 
WRITE(*, 2 
FORMAT(9X, 
WRITE(*,2 
FORMAT(9X, 
WRITE(*,2 
FORMAT(9X, ! 
WRITE(*, 221 
1... 


E а 
WRITE(*,223 
FORMAT(9X, 
WRITE(*,224 
FORMAT(9X, 
WRITE(*, 225 
FORMAT(9X, 
WRITE(*, 226 
FORMAT(9X, 
WRITE(*,227 
FORMAT(9X, 
WRITE(*, 228 
FORMAT(9X, 
WRITE(*,229 
FORMAT(9X, 
WRITE(*,230 
FORMAT(9X, 
WRITE(*,231 
FORMAT(9X, 
WRITE(*,400 
FORMAT(9X,' 
WRITE(*,234 
FORMAT(9X, 
WRITE(*,236 


IMPULSE (TI) WAS TAKEN FROM A 
TIME INTEGRATED ACCELEROMETER. 


TOTAL GAS PULSE IMPULSER TIC) 


WAS OBTAINED FROM FIRING THE 
THRUSTER WITH THE ELECTRODES 
SHORTED. THE DURATION OF THE 


THRUST PULSE (DI) CAME ЕГЕ 


MEASURING THE TIME BETWEEN THE 


HALF CURRENT POINTS ON EOI 


OF CURRENT VS TIME. TH; PR K 


AND FALL OF THE CURRENT WAS 


FAST ENOUGH TO JUSTIFY THE 


SQUARE PULSE SHAPE ASSUMPTION 


INHERENT IN THIS PROCEDURE. 


THE MASS FLOW RATE(DM) WAS 


DETERMINED BY MEARURING THE 


SHAPE... OF THE GAS PULSE UIT 
GAUGE AND COMBINING 


AN ION 


THIS WITH _ THE MASS LOST TRON 
A PLENUM TO THE GA SALT 
CALCULATED DATA WAS OBTAINED 
FROM THE MEASURED DATA USING 
THE FOLLOWING EQUATIONS: 
THRUST(T) = О СЕТА pü 


SPECIFIC IMPULSE SIT (9S DN) 
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Ge 
|") 
|') 
|% 
|") 
|) 


236 FORMAT(9X, EFFICIENCY = T**2/(2*V*J*DM) |") 


! 
| 
WRITE(*,238) 
238 FORMAT(9X,'| AVERAGES OF REPEATED FIRINGS |') 
WRITE(*,239) 
239 FORMAT(9X,'| ARE SENT TO A PLOTTING PROGRAM |') 
WRITE(*,240) 
240 FORMAT(9X,'| WHICH FITS A SPLINED CURVE TO |') 
ШӘШГРЕ(Е 241) XX,YY,ZZ,DMM 
ШІ ЕОЕМАТ(9Х,!| ТНЕМ. ТНЕЦЅТЕВК £',3F4.1,F2.0,' |') 
WRITE(*,242) 
242 FORMAT('+ | 
ПЕ. 0) 


READ(*,243) IK 
243 FORMAT(A1) 
836 CONTINUE 
TET 
L=0 
009 CONTINUE 
WRITE(*,010) I 
010 FORMAT(1X,'ENTER VOLTAGE DATA POINT NUMBER ',I2) 
READ(*,011) V(I) 
011 FORMAT(F5.1) 
WRITE(*,012) I 
012 FORMAT(1X,'ENTER CURRENT DATA POINT NUMBER ',I2) 
READ(*,013) U(I) 
013 FORMAT(F6.0) 
WRITE(*,014) I 
014 FORMAT(IX,'ENTER TOTAL IMPULSE DATA POINT NUMBER ',I2) 
READ(*,015) W(I) 
015 FORMAT(F5.3) 
WRITE(*,950) 

950 FORMAT(1X,'TO CORRECT ANY OF THE LAST THREE VALUES, 
1,'ENTER A 1 AND ENTER THE NEW NUMBERS. OTHERWISE 
2,'ENTER A 0.') 

READ(*,951) JM 

951 FORMAT(11) 
IF(JM.EQ.1) GO TO 009 
WRITE(*,016) 

016 FORMAT(1X,'IF THIS IS THE LAST OR ONLY ELEMENT OF A ' 
1,'GROUPING OF DATA POINTS TO BE AVERAGED TO MAKE ' 
2,'A SINGLE POINT FOR PLOTTING, ENTER THE NUMBER OF' 
3,' ELEMENTS IN THE GROUP. IF IT IS NOT THE LAST' 

4,' ELEMENT IN A GROUP ENTER A 0.') 
READ(*,017) NN 
017 FORMAT(I1) 
IF(NN.EQ.0) GO TO 20 


! 


! 


L=L+1 

VA(L)=0.0 
UA(L)=0.0 
WA(L)=0.0 
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M=I+1-NN 
DO GUS ЛЕШ 
VA(L)=VA(L)+V(J) 
UA(L)=UA(L)+U(J) 
WA(L)=WA(L)+W(J) 

018 CONTINUE 
VA(L)=VACL)/NN 
UA(L)=UA(L)/NN 
WA(L)=WA(L)/NN 
WRITE(*,600) 

600 FORMAT(1X,'THE AVERAGES FOR THIS GROUP ARE') 
WRITE(*,450) VA(L),UA(L),WA(L) 

450 ЕОЕМАТ(1Х,Е10,5,"Ү01Т5!,Р10,2. "АМР5 КТО ЕЕ 
IS ECONDS 

020 CONTINUE 

WRITE(*,30) 
030 FORMAT(1X,'ENTER A 1 IF YOU HAVE MORE DATA POINTS.' 
15 ENTERA ZERO Tr YOUPAPERDOKNE Э) 
READ(*,040) K 
040 ЕОЕМАТ(11) 
IF(K.EQ.0) GO TO (045 
Т=[+1 
СО ТО 009 

045 CONTINUE 
MM=I 
WRITE(*,046) MM,L 

046 FORMAT(1X,'YOU HAVE ENTERED ',I2,' DATA POINTS' 

1,' CONSISTING OF ' M2 CPUS MES 
831 CONTINUE 
WRITE(*,050) 
050 FORMAT(1X,'ENTER THE TOTAL IMPULSE FOR THE GAS PULSE') 
READ(*,060) TIG 
O60 FORMAT(F5.4) 
WRITEC* 070) 
070 FORMAT(1X,'ENTER THE TIME BETWEEN THE ONE-HALF ' 
1, ‘CURRENT POINTS FOR A TYPICAL FIRING OF ' 
2, 'THETHRUSTER') 
READ(*,080) DT 
080 FORMAT(F6.4) 
WRITE(*, 085) 

085 FORMAT(1X,'ENTER A CONTROL-P, THEN ENTER A BLANK TO' 
1,' PRINT THE DATA. AFTER THE DATA IS PRINTED, DO THE' 
2,' SAME THING TO TURN THE PRINTER СЕР 0 

READ(*,089) NP 

089 FORMAT(AI) 

WRITE(*,103) 

103 EORMATL( 26%, 

jon ©) 
WRITE(*“ 104) XX YY 22 


34 


104 FORMAT(45X,'|CATHODE=',F4.1,' CM |ANODE=' 
ee cee | BENCHMARK SCALE=",F4,1,*]”) 
WRITE(*,105) 
105 FORMAT('+ 
! 
B ii 
WRITE(*,350) 
350 FORMAT(46X,' 
I’ 
WRITE(*,351) 
351 FORMAT(45X,'| GAS IMPULSE |THRUST PULSE DURATION'- 
1,'|MASS FLOW RATE |') 
DIETE 352) 
352 FORMAT('+ 
1 
X. ') 
WRITE(*,353) TIG,DT,DM 
553 FORMAT(45X,'|',1X,F5.4,'(N-S)',2X, ' Ох КОН 
Ж 5 [| *,1X,E6.4, 'CKG/S)',2X,'1") 
WRITE(*,354) 
354 FORMAT('+ 
! 
ae ") 
WRITE(*,300) 
300 FORMAT(46X,' 
ЖЕ) 
WRITE(*,087) 
087 FORMAT(45X '| VOLTAGE] AMPERES|TOTAL IMP.(N-S)|' 
ШЕРТЕР | EFRPETENCY(Z) |") 


WRITE(*,500) 

500 FORMAT('+ : 
er: ; 
2 ! > 


D0=2090 I=1,MM 

F(I)=(W(I)-TIG)/DT 

Das 700. 0=E(7)2#72.0/V I )KUCT)*DM*2.0) 
ӨШТІ) -ЕСІ)/(0М%9,8) 

Si 7088) V(T),U(T),W(I),SI(I),E(I) 


ШИНЕКЦИПКАОХ, "| ',ES.1,1X,']',1X,E6.0,'|',4X,F5.3,6X, 
ЖИЕС 0, || ,5Х,ғ4.1,4Х, ||!) 

090 CONTINUE 
WRITE(*,301) 


301 FORMAT('+ 
I 
2 -- ) 
READ(*,093) NA 
093 FORMAT(A1) 
DO 100 IN=1,L 
FA(IN)=(WA(IN)-TIG)/DT 
EA(IN)-100.0*FA(IN)**2.0/(VA(IN)*UA(IN)*DM*2.0) 
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SIA(IN)=FA(IN)/(DM*9.8) 
100 CONTINUE 
WRITE(*,120) 
120 FORMAT(1X,'THE AVERAGE GROUP VALUES ARE.') 
DO 150 I=1,L 
МЕІТЕ(%,130) УА(1), ,UACI), Du ,FACI), DM EA(I) 
130 FORMAT(1X,F5.1,1X, ,1X,F6.0,'|' DADA EB 
Шел SAX EA. a 
150 CONTINUE 
READ(*,700) NB 
700 FORMAT(AI) 
DO 979 III=1,L 
UA(III)=UA(III)/1000. 
979 CONTINUE 
WRITE(*,832) 
832 FORMAT(1X,'ENTER A 1 IF YOU NEED TO ENTER A NEW') 
WRITE(*,833) 
833 FORMAT(1X,'TIME OR MASS FLOW RATE ') 
WRITE(*,834) 
834 FORMAT(1X,'ENTER A 2 TO START OVER, A 3 TO GO ON') 
READ(*,835) NNN 
835 FORMAT(II) 
IF(NNN.EQ.1) GO TO 831 
IF(NNN.EQ.2) GO TO 836 
980 CONTINUE 
WRITE(*,841) 
841 FORMAT(IX,'ENTER 990099 FOR CRT DISPLAY!) 
WRITE(*,842) 
842 FORMAT(1X,'ENTER 965030 FOR A COM2 PLOTTER') 
WRITE(*,985) 
985 FORMAT(1X,'ENTER 960030 FOR A COM1 PLOTTER') 
READ(*,843) IO,MO 
843 FORMAT(I4,I2) 
IF(IO.EQ.9900) 10=10/100 
WRITE(*,972) 
972 FORMAT(IX,'ENTER THE LINE SEGMENT LENGTH!) 
READ(*,973) 2 
973 FORMAT(F6.4) 
VA(L+1)=0.0 
VA(L+2)=60.0 
UA(L+1)=0.0 
UA(L+2)=6.0 
EA(L+1)=0.0 
EA(L+2)=10.0 
SIA(L+1)=0.0 
SIA(L+2)=1000.0 
CALL PLOTS(0,I0,MO) 
CALL WENDOW COs, Oc 1057-725) 
САГЕ РГО оаа) 
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ДР ит“ __— = - _ ___-, „. == ~ ا ا ا‎ а ЕКЕ Br, PASS 


986 
987 


283 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


SYMBOL(-.4,4.3,.16,CHAR(0),90.0,-1) 
SYMBOL(-1.15,4.2,.16,CHAR(2),90.0,-1) 
SYMBOL(-1.9,5.1,.16, CHAR(1), 90.0,-1) 
or O ' = * 90.0,3) 
DINECUA,VA.L, 1991.0) 

CURVE(UA,VA,L,Z) 

LINE(UA,EA,L,1,-1,2) 

CURVE(UA,EA,L,Z) 


LINE( UA, SIA D ] 151) 


GALL CURVE(UA,SIA 22) 

BEE STAXISC.1,.16:1,.05,-1) 

GP PAX IS(0- 0  GURRBNTI(KA)',-13,7.0,0.0,UA(L+1),U 
1A(L+2)) 

CALL AXITS(0.,0.,'VOLTAGE(VOLTS) = ',17,6.0,90.0,VA( 
1L+1),VACL+2)) 

ШЕ ТЕСЕ 75.0. EREEIGIENCY(Z) = ',16,6.0,90.0,EA 


l(L+1),EA(L+2)) 


CALL 


1L+1), 


pau > O N SPECIFIC IMPULSE' 16,6.0 90.0;,SIA( 
SIAC L+2) ) 


ОТА(Ш+2)=500. 
SI(MM+1)=0.0 
SI(MM+2)=500. 
E(MM+1)=0.0 
E(MM+2)=10.0 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


PLOT(0.0,0.0,-999) 

WENDOW(O O. LOT 7.15): 

БОРО СЕ 3 124. =3) 

SYMBOL( =. 443 6 CHAR(2),90.0,-1) 
EINE(ST,E,MM, 1,172) 

CURN HGS Tay oA, by 2) 

ЛЛ 210,.:16,21,..05,-1) 

Ж СЕЕСПЕПЕ IMPULSE”, -16,8.,0,,0.,500.) 
ТШИ НЕ НЕ ПЕНЕР(%2) = 7 165,6.,90:,0.,10.) 
ELOT 002020207999) 


WRITE(*,986) 
FORMAT(1X,'ENTER O IF DONE, 1 IF NOT!) 
READ(*,987) LLL 

FORMAT(II1) 

ШЕШІП 820.0) Go TO 983 

GO TO 984 

CONTINUE 


STOP 
END 
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APPENDIX B 
THEORETICAL MODEL PROGRAM 


PROGRAM MPDI 
REAL J 
F(X)=1.00E-7*(J**4.0/X**2.0-4.94E-12*J**6.0*(X**2.0-2 
1.0*XMIN*X+XMIN**2.0))**0.5 
FF(X)=1.00E-7*(J**4.0/X**2.0-4.94E-12*J**6.0*(XMAX**2. 
10-2.0*XMAX*X+X**2.0))**0.5 

001 CONTINUE 

WRITE(*,010) 

010 FORMAT(1X,'THIS PROGRAM CALCULATES THE SUM OF THE' 
1,' BLOWING, PUMPING, AND ELECTROTHERMAL THRUST 
2,'COMPONENTS FOR AN MPD THRUSTER') 

WRITE(*,005) 
005 FORMAT(1X,'ENTER ANODE LENGTH IN F4.1 CM') 
READ(*,006) Z 
006 FORMAT(F5.2) 
WRITE(*,007) 
007 FORMAT(1X,'ENTER THE MASS FLOW RATE IN KG/S IN F7.5') 
READ(*,008) W 
008 FORMAT(F8.6) 
WRITE(*,020) 
020 FORMAT(1X,'ENTER THE CATHODE RADIUS IN F4.1 CM') 
READ(*,030) XMIN 
030 FORMAT(F5.2) 
WRITE(*,040) 
040 FORMAT(1X,'ENTER THE ANODE RADIUS IN F4.1 CM') 
READ(*,050) XMAX 
050 FORMAT(F5.2) 
056 FORMAT(F7.1) 
WRITE(*,060) 
060 FORMAT(1X,'ENTER AN EVEN NUMBER OF INTERVALS IN I100') 
READ(*,070) N 
N-10 
061 CONTINUE 
WRITE(*,055) 
055 FORMAT(1X,'ENTER AMPERES IN F7.1') 
READ(*,056) J 
070 FORMAT(I100) 
H=(XMAX-XMIN)/(2.0*N) 
A=0.0 
AA=0.0 
AAA=0.0 
X=XMIN 
XXX=(XMAX+XMIN)/2.0 
DO 100 I=1,N 
ШЕТІН 
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131 
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133 


200 


134 


201 


222 


140 


150 


1521 


XXXX=XXX+H 
A=A+H*(F(X)+F(XX))/2.0 

AA=AA+H*(FF(XXX)+FF(XXXX))/2.0 

X=XX 

XXX ос 

CONTINUE 

AAA=A+AA 

WRITE(*,130) AAA 

FORMAT(1X,'BLOWING THRUST = ',F8.4) 

A=XMAX 

B=XMIN 

Y=(XMIN+XMAX)/2.0 
AAAA=4.44E-13*(7,0*Y**3,0/6.0-(B**3.,0+A**3.0)/3.0-(A*B 
1**2.04B*A**2.0)/4.0-(Y**2.0/2.0) * (A*FLOG(Y/A) xB*LOG(Y/B 
2) )+(A*B**2 .O*LOG(B/A)+B*A**2.0*LOG(A/B))/2.0)*J**3.0/Z 
WRITE(*,131) AAAA 7 | 


FORMAT(1X,'PUMPING THRUST =',F8.4) 
AAAAA=3.25E-7*J**1,5*W**0,3*3,14159*(XMAX**2.0-XMIN**2 
1720) 

WRITE(*,132) AAAAA 

FORMAT(1X,'ELECTROTHERMAL THRUST =',F8.4) 
AX=AAA+AAAA+AAAAA | 

WRITE(*,133) AX | 

FORMAT(1X,'TOTAL THRUST =',F8.4) 


Al=(55.4-26.31*#( I-XMAX/5. De 22.27 7-77 Dun 
B1=(39.214-65.84*(1-XMAX/5:1)-27.656*(1-Z/21.6))*1E-4 
C1=(35.7-54.5*(1-XMAX/5.1)+3.88*(1-Z/21.6))*1E-8 
V=(A1-B1*J+C1*J**2)*(178.-47.3*1000. *W+4.3*( 1000. *W) 
1**2)/49. 

WRITE(*,200) V,J 

FORMAT(1X,'VOLTAGE =',F9.4,2X,'CURRENT =',F7.1) 
YX=AX/(W*9.8) | 

WRITE(*,134) YX 

FORMAT(1X,'SPECIFIC IMPULSE =',F8.3) 
E=(.5*W*(YX#9,8)**2.0)/(V*J) 

WRITE(*,201) E 

FORMAT(IX,'EFFICIENCY =',F8.4) 

P=V¥J 

WRITE(*,222) P 

FORMAT(1X,'POWER =',F9.1) 

WRITE(*,140) 

FORMAT(1X,'ENTER 1 TO START OVER, 2 TO CHANGE CURREN' 
iT ONLY S SISA 

READ(*,150) M 

FORMAT(I1) 

IF(M. EQ. 1) CO TO 001 

ПЕСИ ЕО CO TO 061 

CONTINUE 


STOR 
END 
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